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We discuss the role of the long-range elastic interaction between the contacts inside an inho- 
mogeneous frictional interface. The interaction produces a characteristic elastic correlation length 
A c = a? Ejk c (where a is the distance between the contacts, k c is the elastic constant of a contact, 
and E is the Young modulus of the sliding body), below which the slider may be considered as a 
rigid body. The strong inter-contact interaction leads to a narrowing of the effective threshold distri- 
bution for contact breaking and enhances the chances for an elastic instability to appear. Above the 
correlation length, r > A c , the interaction leads to screening of local perturbations in the interface, 
or to appearance of collective modes — frictional cracks propagating as solitary waves. 
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I. INTRODUCTION 

Studies of sliding friction, a subject with great prac- 
tical importance and with rich physics, attracted an in- 
creased interest during last two decades [HQ. Tip-based 
experimental techniques as well as atomistic molecular 
dynamics (MD) computer simulations describe with con- 
siderable success the processes and mechanisms operat- 
ing in atomic-scale friction. Much less is known when 
going to meso- and macro-scale friction, where one has 
to take into account that the frictional interface is inho- 
mogeneous and generally complex. An immediate exam- 
ple is dry friction between rough surfaces. Even when 
the sliding surfaces are ideally flat but for example the 
substrates are not monocrystalline, or there is an inter- 
posed solid lubricant film consisting of misoriented do- 
mains, the frictional interface is again inhomogeneous. 
The same may be true even for liquid lubrication, if un- 
der applied load the lubricant solidifies making bridges 
due to Lifshitz-Slozov coalescence. In these cases, the 
so-called earthquake-like (EQ) type models can be suc- 
cessfully applied [2h13||. In the EQ model, the two (top 
and bottom) mutually sliding surfaces are coupled by a 
set of contacts, representing, e.g., asperities, patches of 
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lubricant, or 2D crystalline domains. A contact is as- 
sumed to behave as a spring of elastic constant k c so long 
as its length is shorter than a critical value x s — f s /k c ; 
above this length the contact breaks, to be subsequently 
restored with lower stress. The sliding kinetics of this 
model may be reduced to a master equation (ME), which 
allows an analytical study @, [TJ [lj| . 

In the simplest approach, the slider is treated as a rigid 
body. Due to the non-rigidity of the substrates, however, 
several length scales naturally appear in the problem. 
First, different regions of the interface will exhibit differ- 
ent displacements. The length Xl such that for distances 
r ^ Xl the displacements are i ndep endent, is known as 
the Larkin-Ovchinnikov length |l4{ . It was shown [l5[ 
that for the contact of stiff rough solid surfaces Xl may 
reach unphysically large values ~ i0 100000 m Second, 
deformation of the solid substrates leads to the elastic 
interaction between the contacts. Elasticity will corre- 
late variations of forces on nearest contacts over some 
length A c known as the elastic correlation length [Tfjl ]. 
Third, displacements in one region of the slider will be 
felt in other regions on the distance scale set by of a 
screening length X s . Finally, the breaking of one contact 
may stimulate neighboring contacts to break too (the so- 
called concerted, or cascade jumps), following which an 
avalanche-like collective motion of different domains of 
the interface may appear [f|. 

In this paper we discuss collective effects in the fric- 
tional interface and propose approaches to treat them 
from different viewpoints. In particular, our aim is to 
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clarify the following questions: (i) what is the law of in- 
teraction between the contacts; (ii) at which scale can the 
slider be considered as a rigid body, or what is the coher- 
ence distance A c within which the motion of contacts is 
strongly correlated; (Hi) whether the interaction effects 
can be incorporated in the master equation approach; 
(iv) how does the interaction modifies the interface dy- 
namics; (v) what is the screening length A s ; (vi) when 
do avalanche motion of contacts (a self-healing crack) ap- 
pear, and what is the avalanche velocity? 

The paper is organized as follows. The earthquake- 
like model, its description with the ME approach, and the 
elastic instability responsible for the stick-slip motion are 
introduced in Sec. II. The interaction between contacts 
is studied in Sec. IIIII An approach to incorporate the 
interaction between contacts into the ME approach in 
a mean-field fashion is described in Sec. IIV1 The role 
of interaction at the meso/macro-scale is considered in 
Sec. El Finally, discussions in Sec. [VI] conclude the paper. 



II. EARTHQUAKE-LIKE MODEL, MASTER 
EQUATION AND ELASTIC INSTABILITY 

A. The earthquake-like model 

In the EQ model the sliding interface is treated as a set 
of N contacts which deform elastically with the average 
rigidity k c . The zth contact connects the slider and the 
substrate through a spring of shear elastic constant hi. 
When the slider is moved, the position of each contact 
point changes, the contact spring elongates (or shortens) 
so that the slider experiences a force —F — ^ /i from 
the interface, where fi = kiXi and Xi(t) is the shift of the 
ith junction from its unstressed position. The contacts 
are assumed to be coupled "frictionally" to the slider. 
As long as the force \fi\ is below a certain threshold f S i, 
the zth contact moves together with the slider. When 
the force exceeds the threshold, the contact breaks and a 
rapid local slip takes place, during which the local stress 
drops. Subsequently the junction is pinned again in a 
less-stressed state with fbi , and the whole process repeats 
itself. Thus, with every contact we associate the thresh- 
old value f s i and the backward value , which take ran- 
dom values from the distributions P c (f) and R(f) corre- 
spondingly. When a contact is formed again (re-attached 
to the slider), new values for its parameters are assigned. 
The EQ model was studied numerically in a number of 
works [3-[l0|. typically with the help of the cellular au- 
tomaton numerical algorithm. 



B. The master equation approach 

Rather than studying the evolution of the EQ model 
by numerical simulation, it is possible to describe it an- 
alytically [1, [HI, [Hj] . Let P c (x) be the normalized prob- 
ability distribution of values of the stretching thresh- 



olds x S i at which contacts break; it is coupled with 
the distribution of threshold forces by the relationship 
P c (x) dx — P c (f) df, i.e., the corresponding distributions 
are coupled by the relationship P c (x) oc xP c [f(x)], where 
/ oc x 2 [ll|- We assume that the distribution P c (x) has 
a dispersion Ax s centered at x — x c . 

To describe the evolution of the model, we introduce 
the distribution Q(x;X) of the contact stretchings Xi 
when the sliding block is at position X. Evolution of the 
system is described by the integro-differential equation 
(known as the master equation, or the kinetic equation, 
or the Boltzmann equation) [H, [ll| 
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dX dx 



+ P{x) 



Q(x;X) = R(x)T(X), (1) 



where 



/oo 
d£P{QQ(t;X) 
-oo 



(2) 



and 



P(x) = P c (x)/J c (x) , J c (x) 



d£Pc(Z). (3) 



Then, the friction force (the total force experiences by 
the slider from the interface) is given by (k c = (hi)) 



/oo 
dxxQ(x;X) 
-OO 



(4) 



In the steady state corresponding to smooth sliding, 
the ME reduces to 



dQ{x)/dx + P(x) Q{x) = R(x) T , 
which has the solution 



,{x)=ME P (x) 



1 + T f d£R(£)/E P (t) 
Jo+ 



(5) 



(6) 



where J\f is the normalization constant, dxQ s (x) = 1 , 
and 



E P {x) = exp {-U(x)} , U(x) = f (%P(Z) . 

Jo 



(7) 



C. Elastic instability 



The solution of the ME [H, [TTJ] shows that when a rigid 
slider begins to move adiabatically, X > 0, it experiences 
from the interface a friction force F co (X) < 0. Initially 
|Foo| grows roughly linearly with X, \Foo \ ss K S X (here 
K s = Nk c is the total elastic constant ( "rigidity" ) of the 
interface), until it reaches a value ~ F s — AF S , where 
F s « K s x c and A_F S « K s Ax s . Gradually however con- 
tacts begin to break and reform, slowing down the in- 
crease of | -Fool and then inverting the slope through a 
displacement Ax s until almost all contacts have been 
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reborn. Successively the process repeats itself with a 
smaller amplitude until, due to increasing dispersion of 
breaking and reforming processes, the force asymptoti- 
cally levels off and attains a position independent steady 
state kinetic friction value with smooth sliding. 

According to Newton's third law, the external driving 
force Fd — K(vt — X) which causes the displacement X 
(here K is the slider rigidity and v is the driving velocity), 
is compensated by the force from the interface, Fd = 
F(X). Smooth sliding is always attained with a rigid 
slider. It persists for a nonrigid slider as well, so long 
as the pulling spring stiffness is large enough, K > K* , 
where 

K* = maxi^pO , F^(X) = dF OQ (X)/dX . (8) 

When conversely the slider, or the pulling spring elastic 
constant are soft enough (K < K*) there is a mechanical 
instability. The driving force Fd cannot be compensated 
by the force from the interface, and the slider motion be- 
comes unstable at X c , where X c is the (lowest) solution 
of F^(X) = K (for details see Refs. [| EH). The me- 
chanical instability yields stick-slip frictional motion of 
the slider. 

Thus, the regime of motion — either stick-slip for K <C 
K* or smooth sliding for K 3> K* — is controlled by the 
effective stiffness parameter K* ~ K s x c / Ax s . When all 
contacts are identical, Ax s — so that K* = oo, one 
always obtains a stick-slip motion. 



p = 10 g/cm 3 and <r c = 10 9 N/m 2 (steel), we obtain 
72 ~ 10~ 3 which should be typical for a contact of rough 
stiff surfaces. For softer materials, and especially for a lu- 
bricated interface, the values of 72 would be much larger, 
e.g., 72 ~ 0.1. 

The second dimensionless parameter 71 = k c /Ea char- 
acterizes the stiffness of the contacts. To estimate 71, 
assume again contacts with the shape of a cylinder of 
radius r c and length h (h is the thickness of the in- 
terface). Suppose in addition that one end of a con- 
tact ("column") is fixed, and a shear force / is ap- 
plied to the free end. This force will lead to the dis- 
placement x — f/k c of the end, where k c = 3E c I/h 3 , 
E c is the Young modulus of the contact material and 
I = Trr^/A is the moment of inertia of the cylinder [17]. 
In this way we obtain k c = (3tt / 1 A)(E c r c )(r c /h) 3 , so 
that 71 = (3tt/4) (E c a 3 /Eh 3 ) {r c jaf = 707! with 7o = 
(3ir 1 4)(E C / E)(a/ h) 3 . For the contact of rough surfaces, 
where E c — E and a > h, we have 70 ^ 1, while for 
lubricated interfaces where E c <C E, one would expect 
7o ^ 1- 

An estimate of characteristic values [l[ leads to r c ~ 
(10 -3 -7- 10~ 2 ) a. Thus, for the steel slider considered 
above, taking r c = h = 1 /im and intercontact spacing 
a = 3 x 10 2 r c , we obtain N ~ 10 3 and k c ~ 5 x 10 5 N/m, 
so that the global stiffness of the interface is K s ~ 5 x 
10 s N/m. 



D. Material parameters 

It is useful here, before proceeding with the analyti- 
cal and numerical developments necessary to answer the 
questions posed in the Introduction, to review the practi- 
cal significance and magnitude of the model parameters. 

Elastic constant of the slider. The slider (shear) elas- 
tic constant K is equal to K = [E/2 (1 + a)][L x L y /H], 
where L x , L y and H are the slider dimensions, E and a 
are the substrate Young modulus and Poisson ratio, re- 
spectively 17J. For example, for a steel slider of Young's 
modulus E — 2 x 10 11 N/m 2 , Poisson's ratio a = 0.3 and 
the size L x x L y x H — 1 cm x 1 cm x 1 cm, we obtain 
K - 10 9 N/m. 

Rigidity of the interface contacts. Here we characterize 
the typical magnitudes of the contact stretching length 
x c and stiffness k c . Assume the slider and the substrate 
to be coupled by N = L x L y /a 2 contacts, and that the 
contacts have a cylindrical shape of (average) radius r c 
with a distance a between the contacts. It is useful to 
introduce the dimensionless parameter 72 = r c /a, which 
may be estimated as follows (l| . Consider a cube of linear 
size L on a table. The weight of the cube Fi — pL 3 g (p 
is the mass density and g = 9.8 m/s 2 ) must be compen- 
sated by forces from the contacts, Fi = Nr^a c , where a c 
is the plastic yield stress. Then, 7I = (Nr^)/(Na 2 ) — 
(pL 3 g)/(a c L 2 ), or 72 = (pLg/a.) 1 / 2 . Taking L = 1 cm, 



Stick-slip versus smooth sliding. As mentioned above 
in Sec. Ill CI the regime of motion (either stick-slip or 
smooth sliding) is controlled by the parameter K* ~ 
K s x c / Ax s . For the steel slider considered above, esti- 
mates gave K - 10 9 N/m and K 8 ~ 5 X 10 8 N/m. Thus, 
if the surfaces are rough so that Ax s ~ x c , then K > K* 
and one should typically get smooth sliding. Stick slip 
appears further disfavored if we consider a realistic P c {x) 
distribution. For all cases mentioned in Introduction — 
the contact of rough surfaces (both for elastic or plastic 
asperities), the contact of polycrystal (flat) substrates, 
and the case of lubricated interface, when the lubricant, 
melted during a slip, solidifies and forms bridges at stick, 
— the distribution P c [x) is rather wide with a large con- 
centration of small-threshold contacts ll| , which makes 
the value of K* very small. Thus, the theory predicts 
that most systems do not undergo an elastic instability 
and should not therefore exhibit stick-slip. This conclu- 
sion contradicts everyday experience as well as careful 
experiments, where stick-slip is pervasive. As suggested 
by EQ simulations [13j |. the discrepancy is most likely 
caused by ignoring the elastic interaction between the 
contacts. 

The role of interaction is considered in the next sec- 
tions. First, however, we need to define the form and 
parameters of the interaction between contacts. 
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FIG. 1. Left: decaying of the displacement field at the inter- 
face (schematic): (a) for a single contact u(r) oc r _1 , (b) for 
a single hole u(r) oc r~ 3 , and (c) for the array of contacts. 
Right: change of forces on contacts when the central contact 
is removed (71 = 0.06). 



III. INTERACTION BETWEEN CONTACTS 

Friction is not a simple sum of individual contact prop- 
erties. The collective behavior of the contacts is im- 
portant. Recently Persson fl8l - l22j developed a contact 
mechanics theory based on the fractal structure of sur- 
faces in order to determine the actual contact area at 
all length scales, which determines the friction coeffi- 
cient. This approach includes the presence of multi- 
ple contacts and leads to the correct low-threshold limit 
P c (f —> 0) = 0. Persson found that the distribution 
of normal stresses a (a > 0) at the interface may ap- 
proximately be described by the expression P a {cr) oc 
exp [— (a — a) 2 /Aa 2 ] — exp [-(a + a) 2 /Aa 2 ] , where a 
is the nominal squeezing pressure, the distribution width 
is given by Act = E*TZ 1/2 (E* is the combined Young 
modulus of the substrates, E*^ 1 = E^+E^ 1 where E\^ 
are the Young modula of the two substrates) , and the pa- 
rameter 1Z is determined by the roughness of the contact- 
ing surfaces, K = ^tt)" 1 / dq q 3 f d 2 x (/i(x)/i(0))e- iqx . 
Assuming that a local shear threshold is directly pro- 
portional to the local normal stress, / cx a, we finally 
obtain the distribution, which is characterized by a low 
concentration of small shear thresholds, P c (f) oc f at 
/ 0, and a fast decaying tail, P c (f) oc exp(— f 2 /f* 2 ) 
at / — > 00, i.e., now the peaked structure of the distribu- 
tion is much more pronounced. 

However an important aspect which has to be included 
is the redistribution of the forces when some contacts 
deform or break. A concerted motion of contacts may 
emerge only due to interaction between the contacts 
which occurs through the deformation of the bulk in the 
directions parallel to the average contact plane. It is this 
aspect that we want to consider here. For the elastic 
interaction, a qualitative picture is presented in Fig. Q] 
(left). When a contact acts on the surface at r = with 
a force /, it produces a displacement field u(r) oc 
which affects other contacts (Fig. [TJl) — similar to the 



Coulomb potential for a point charge [17|. However, if 
there are two surfaces, then the same contact acts on 
the second surface with the opposite force — / and, if the 
two surfaces are in contact, the resulting displacement 
field should fall as u(r) oc r~ 3 (Fig. [T|d) — similar to the 
dipole-dipole potential for a screened point charge near a 
metal surface [2j|. The question thus is the form of the 
interaction for the multi-contact interface (Fig. [TJ;). We 
will show that the interaction between the contacts has a 
crossover from the r _1 slow Coulomb decay at short dis- 
tances to the faster dipole-dipole one at large distances. 

A. Analytics 

Let us consider an array of N elastic contacts (springs) 
with coordinates = {xi,yi,0}, i — 1,...,N, between 
the two (top and bottom) substrates. If the interface 
is in a stressed state, the contacts act on the top sub- 
strate with forces ii = {fix, fiy, fiz}- These contact 
forces produce displacements u^ top ' > of the (bottom) sur- 
face of the top substrate. The 3-/V-dimensional vectors 
U(top) = {Uj-* 015 -'} and F t = {£;} are coupled by the linear 
relationship Ij( top ) = G( top )F t . Elements of the elastic 
matrix G' top - ) (known also as the elastic Green tensor) for 
a semi-infinite isotropic substrate were given by Landau 
and Lifshitz [l7| : 

G ix , jx = g(n 3 )[2(l - a) + 2ax%lr%\ 

Gix,jy = 2g( r ij) <Jx ijVijl r ij 

Gix, ]Z = -g(nj)(i - 2a) Xij/nj (9) 

Giz, ]Z = 2flr(ry)(l - a) , 

where Xij — Xi — Xj, g(r) = (1 + a) / {2n Er) , and a and 
E are the Poisson ratio and Young modulus of the top 
substrate, respectively. 

In the equilibrium state, the forces that act from the 
contacts on the bottom substrate, must be equal to 
Fb = — F t according to Newton's third law. These forces 
lead to displacements of the (top) surface of the bottom 
substrate, Tj( bottom ) = _G( bottom )Ft. The elements of 
the bottom Green tensor G( bottom ) are defined by the 
same expressions © except the xz elements for which 

G-k° t * om * ) = — G^jl (if the substrates are identical, the 
z displacements are irrelevant). Thus, the relative dis- 
placements at the interface due to elastic interaction be- 
tween the contacts are determined by 

U _ uCtop) _ uOottom) = _q f ) ( 10 ) 

where F = -F t and G = G< top ) + G( bottom ). 

On the other hand, the forces and displacements are 
coupled by the diagonal matrix (the contacts' elastic ma- 
trix) K, Kiajp = k ia 6 i: j5 al 3 (a, (3 = x,y,z): 

F = K(U + U), (11) 

where Uo defines a given stressed state (because of linear- 
ity of the elastic response, final results should not depend 
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of Uo). The total force at the interface, f = J2i fii must 
be compensated by external forces applied to the sub- 
strates, e.g., by the force f( cxt ) = f applied to the top 
surface of the top substrate if the bottom surface of the 
bottom substrate is fixed. 

Combining Eqs. (JTOj) and ([TT]l. we obtain F = K (U - 
GF), or 

F = BKU , where B = (1 + KG)" 1 . (12) 

If one changes the contact elastic matrix, K — >• K+5K, 
then the interface forces should change as well, F — > 
F + SF. From Eq. we have SF = (<5B)KU + 

B((JK)Uo. Then, <5B may be found from the equation 
S[B (1 + KG)] = (<5B)(1 + KG) + B(£K)G = 0. There- 
fore, finally we obtain: 

<5F = B£K(1-GBK)U . (13) 

Above we have assumed that <5K is small. If it is not 
small, we have to use the expression SF = B SK (1 — 
GBK)U , where K = (1 + <SKGB) 1 (K + SK). 

Now, if we remove the i*th contact by putting 5ki a — 
—kiaSii* and then calculate the resulting change of forces 
on other contacts, we can find a response of the interface 
to the breaking of a single contact as a function of the 
distance r = — r> from the broken contact. 

B. Numerics 

Equation (|13p may be solved numerically by standard 
methods of matrix algebra. We explore an idealized ar- 
ray of identical contacts, ki a = k c and (Vo)ia — uoS ax 
for all i, organized in a square 89 x 89 lattice with spac- 
ing a = 1, with the broken contact i* at the center of the 
lattice. For singular terms of the Green function © we 
apply a cutoff at ru = r c . Numerical results depend on 
two dimensionless parameters. The first is 71 = k c / 'E*a, 
which determines the stiffness of the array of contacts rel- 
ative the substrates (here E^ 1 = E^ + E^ ttom ). The 
second parameter 72 = r c /a characterizes a single con- 
tact (or the density of asperities). For the Poisson ratio 
we took a typical value a — 0.3. A typical distribution of 
breaking induced force changes is shown in Fig. Q] (right). 

The numerical results for the x-component of dimen- 
sionless force Sf = 5F x /(k c u Q ) are presented in Fig. [2] 
The function Sf(r) exhibits a crossover from a slow 
Coulomb like decay 5f(r) cx r~ x at short distances r«A c 
to the fast dipole-dipole like decay Sf(r) cx r~ 3 at large 
distances r ^> A c . The near and far zones are sep- 
arates by the elastic correlation length A c first intro- 
duced by Caroli and Nozieres [16|. It may be estimated 
in the following way: the stiffness of the "rigid block" 
K ~ EX C should be compensated by that of the inter- 
face, K ~ k c (A c /a) 2 (stiffness of one contact times the 
number of contacts). This leads to 

A c « a/71 = a 2 E/k c . (14) 




X 

FIG. 2. (Color online): Dependence of the change of forces 
Sf(r) on the distance x from the broken contact for three val- 
ues of the interface stiffness: 71 = 0.003 (blue down triangles, 
dashed line), 0.06 (red solid circles, dotted line) and 0.8 (black 
up triangles, solid line) at fixed value of 72 = 0.3 (a = 0.3). 
The lines show the corresponding power laws. 

The rigid slider corresponds to the limit E — > 00, or 
71 — > 0. Therefore, the slider may be considered as a 
rigid body (e.g., in MD simulation), if its size is smaller 
than A c . For the steel slider considered in Sec. Ill D[ es- 
timation gives A c / a ~ 10 2 . Up to distance A c the con- 
tacts strongly interact. If the ith contact breaks and 
its stretching changes on \Sxi\ ~ x c , then the force on 
the jth contact at a distance ry < A c away, changes by 
Sfj w kk c aSxi/rij, where the dimensionless parameter 
k < I characterizes the strength of interaction (numerics 
gives k ~ 10 -3 ). In the near zone r< A c the interaction 
between the contacts may be accounted for within the 
master equation approach in a mean-field fashion as de- 
scribed in the next Sec. HVl At larger distances, different 
regions of the slider will undergo different displacements. 
Therefore, in the far zone, r > A c , wc must take into 
account the elastic deformation of the slider. 



IV. NEARBY CONTACTS: MEAN FIELD 
APPROACH 

EQ model with interaction between the contacts. Let 
us now include the dynamical interaction between the 
contacts. When a contact breaks, the now unsustained 
shear stress must be redistributed among the neighbor- 
ing contacts. We assume that because of elastic interac- 
tion between the contacts i and j, the forces acting on 
these contacts have to be corrected as fi —> /, — A/y 
and fj -> fj + Afij, where A/ y = fey (xj - x») in lin- 
ear approximation. For example, let at the beginning 
the contacts be relaxed, Xj(0) — #i(0) = 0. Due to slid- 
ing motion, all stretchings grow together, so that still 
A/y = 0. At some instant t let the jth contact break, 
Xjit) — > 0, with the ith contact still stretched, Xi(t) > 0. 
Clearly, as the jth contact breaks, the force on the zth 
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FIG. 3. (color online): The steady state distribution Q s (x) for 
the rectangular threshold distribution P c o(x) with x s = 1 and 
Ax 3 = 0.25 and different values of the interaction strength 
K = 0, 0.02, 0.06, 0.1, and 0.5. The EQ simulations (dotted) 
are compared with the ME results (solid curves). 
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FIG. 4. (color online): Onset of sliding: the initial part of the 
dependence of the friction force F on the slider displacement 
X for different strength of interaction k = 0.005 (black), 0.01 
(cyan), 0.03 (red), 0.05 (blue), and 0.07 (magenta). Dotted 
curves show the results of EQ simulation, and solid curves, the 
mean- field ME approach. The threshold distribution P c o(x) 
has the rectangular shape with x s = 1 and Ax s = 0.25. 



contact increases, A/y(i) = —fey Xi(t) < 0. The ampli- 
tude of interaction decreases with the distance r from the 
broken contact as Af oc r" 1 at short distances r < A c . 
Neglecting the anisotropy of interaction, we assume that 
fey = f/\rij\, where / is a parameter. 

We simulated a triangular lattice of N = 60 x 68 = 
4080 contacts with periodic boundary conditions and lat- 
tice constant a = 1, with an average contact spring con- 
stant k c = 1 and radius of interaction A c = 3a or A c = 5a. 
We assumed fu = and a rectangular shape of the dis- 
tribution P c (x), i.e., P c (x) — P c o(x) — (2A# S ) _1 for 



\x — x s \ < Ax s and otherwise, which admits an exact 
solution for noninteracting contacts [ll| (more realistic 
distributions give the same results). 

Figures[3]and|3]show the result of simulations for differ- 
ent values of the dimensionless strength of the interaction 

k = f/(k c x c ) , (15) 

where x c = J dxxP c o(x) is the average stretching of 
the initial threshold distribution (for the rectangular dis- 
tribution x c = x s ). These results yield the following 
conclusions. First, in the steady state, the interaction 
causes a narrowing of the final distribution Q s (x). At 
high interaction strength k, the distribution approaches 
a narrow Gaussian. Second, the drop of frictional force 
F(X) at the onset of sliding (at X ~ x c ) gets steeper 
and steeper as k grows. Therefore, contact interactions 
reinforce elastic instability. Third, above a critical inter- 
action strength, k > k c ~ 0.1, a multiplicity of contacts 
break simultaneously at the onset of sliding, and there 
is an avalanche, where the force F (X) drops abruptly. 
The average avalanche size may be estimated similarly 
as done in Ref. 0. 

While the full EQ model may be only studied numeri- 
cally, it is always useful to have analytical results, even if 
only of qualitative level. In what follows we show that the 
main EQ results may be reproduced within the ME ap- 
proach by using "effective" P c (x) and R{x) distributions 
defined in a mean-field fashion. In this section the ME 
equation is only used to reproduce the EQ results. This 
is however useful because it provides an additional un- 
derstanding of the results as the effective distributions, 
obtained in this analysis, provide a description of the 
collective effects affecting the contacts in terms of simple 
functions. 

Smooth sliding. Using the steady state solution of the 
ME, Eqs. ([6]) and (O, one may approximately recover 
the functions P c (x) and R{x) if the stationary distribu- 
tion Q s {x) is known. Indeed, for small x, where P(x) 
is close to zero, the left-hand side of Q s (x) allows us to 
find R(x) as R(x) oc Q' s (x) (see Eq. ([5])), while the right- 
hand side of Q s (x), where x ~ x c and the contribution 
of R(x) to the shape of the steady state distribution is 
negligible, gives us [ll| P c (x) oc P(x)Q s (x) oc —Q' s (x). 
Thus, differentiating the function Q s (x) obtained in the 
EQ simulation, we may guess shapes of the effective dis- 
tributions P c {x) and R(x) which, when substituted in 
the ME, would produce a solution Q s (x) close to that 
obtained in the EQ simulation. 

Using the simulation results, let us suppose that the de- 
tached contacts form again with nonzero stretchings, i.e., 
that the distribution R(x) is shifted to positive stretching 
values, 

R(x) = G(x - ax c , 7£ c ) , (16) 

where G(x, a) is the Gaussian distribution with zero 
mean and standard deviation a, 

G(x i a) = -L=e X p(~). (17) 
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At the same time, we suppose that the effective threshold 
distribution P c {x) shrinks and shifts with respect to the 
original ( "noninteracting" ) one, 

P h {x) = (3P c0 \p{x - ax c )\ . (18) 

Let us moreover take its convolution with the Gaus- 
sian function (|T7|). P c (x) = P h ® G = J d£ P h (x - 

The results of this procedure for the rectangular dis- 
tribution P c0 (x) are shown in Fig. [3] We see that with a 
proper choice of the parameters a, (3 and 7, the ME so- 
lutions Q s (x) perfectly fits the numerical EQ results (for 
the parameters a, f3 and 7 in Fig. [3] we used expressions 
P = 1 + b\K, a = 62K//3 and 7 = 63a — 64a 2 with the co- 
efficients h = 18, 6 2 = 9.6, 63 = 0.142 and 64 = 0.232). 
Results of similar quality were also obtained for other 
simulated cases, e.g., for larger radius of the interaction 
or for wider threshold distribution P c q(x). 

The dependences of the fitting parameters a, /3 and 
7 on the dimensionless strength of interaction k may be 
found in the following way. To begin with, for nonin- 
teracting contacts initially a = 7 = and (8 = 1. It 
is reasonable to expect that in the lowest approxima- 
tion a, 7 oc k and /? — 1 cx n. Indeed, because the shift 
of the effective distributionP c (x) appears because of the 
interaction, af c = ■ Afij , at small n we have approx- 
imately 

a ~ 0.5 a -2 / d 2 r Kx c /\r\ — ttkXcXc/o 2 . (19) 
Jo 

At large k, however, a has to saturate, e.g., as a oc k/0, 
because the shift cannot be larger than x c , i.e., a < 1. 
Then, because the distribution P c (x) shrinks from both 
sides, we have b\ ~ 2 & 2 - 

Thus, the interaction makes the threshold distribution 
P c (x) narrower by a factor /3 and shifts its center to the 
left-hand side, x c — > vx c , where v = a + changes 
from 1 to 0.5 as the interaction strength k increases from 
zero to infinity. 

Onset of sliding. The beginning of motion when 
started from the relaxed configuration, Q(x; 0) = 5(x), 
cannot be explained by the approach used above, be- 
cause the effective distribution P c (x) is "self-generated" 
during smooth sliding, i.e., it can be applied only when 
the process of contacts breaking-reattachment is contin- 
uously operating. Nevertheless, the initial part of the 
F(X) dependence may still be described by the effective 
ME approach, but with the modified "forward" threshold 
distribution given by the expression 

P cl {x) = Nx e °P c0 [0 Q (x - a x c )\ , (20) 

where N is a normalization factor, dx P C i(x) = 1. 
The parameter ckq is now defined so as to keep the lowest 
boundary unshiftcd, /?o(^fixO — ctox c ) = Xfi x o witha^o = 
x L = x s - Ax s , so that a = (x &x0 /x c )(l - /3' 1 ). The 
"backward" distribution R(x) is still defined by Eq. (fTBT) 
with the same parameters as above. 




1 E 1 I 1 I 1 I 1 I . I 1 I 1 L_ 

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 



K 

FIG. 5. The effective interface stiffness K* s (normalized on 
the noninteracting value) as a function of the strength of in- 
teraction k for the realistic threshold distribution P c o{x) = 
(2/x s )u 3 e~ u2 , u = x/x 3 with x s = 1, when K*/K 3 = 0.179. 

Numerics shows that with a proper choice of the fit- 
ting parameters /3q and eo for a given value of n, the ini- 
tial part of the function F(X) may be reproduced with 
quite high accuracy. Moreover, for a rather wide range of 
k values, the EQ simulation results may be reproduced 
by the ME approach with a reasonable accuracy using 
only three fitting parameter c\, c 2 and k c , if the param- 
eters (3q and eo in Eq. (|2"0")) are given by the expressions 
/3o = 1 + c\k/(1 — k/k c ) and eo = C2(/3o — 1), where the 
parameter k c corresponds to the critical "breakdown" in- 
teraction strength when many contacts begin to break 
simultaneously. For k > k c , the drop of F(X) becomes 
jump-like, so that K* = 00 and stick-slip will appear for 
any stiffness of the slider K < 00. Note that the value of 
k c may be estimated from the equation ax c ~ Ax s . 

For the rectangular shape of the distribution P c o(x) 
the result of this procedure is demonstrated in Fig. [H 
(the fitting parameters are c\ — 33.9, c 2 = 3.0 and k c — 
0.074). 

Of course, the P C i(x) function, Eq. (|2"0")) , can describe 
only the initial part of the F(X) dependence, when F(X) 
grows, reaches the first maximum and then decreases. To 
simulate the whole dependence F(X), one would have to 
involve the evolution of P c (x) with sliding distance, e.g ., 
as some "aging" process P CJ (x) — > P c (x) (see Ref. \U\ ) 
with the initial distribution P c .im(x) — P c i{x) and the 
final one P c ,fm[%) — Pc(x). 

Stick-slip versus smooth sliding. As was mentioned 
above, stick-slip appears as a result of elastic instabil- 
ity which is controlled by the relation between the slider 
stiffness K and the effective interface stiffness K* . For 
noninteracting contacts K* w K s x c / Ax s ; because typ- 
ically Ax s ~ x c , estimates give K* < K so that stick- 
slip should never appear. The interaction between con- 
tacts strongly enhances the elastic instability thus mak- 
ing stick-slip much more probable. Indeed, because of 
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the effective shrinking of the threshold distribution, the 
parameter K* increases roughly as K* — > K* s ~ PqK* , 
i.e., the effective interface stiffness K* s grows with the 
strength of interaction n, and the elastic instability can 
now appear. For example, for a realistic threshold distri- 
bution the dependence of K* s on the strength of inter- 
action k is shown in Fig. [SJ 

The strength of interaction between the contacts may 
be found as k « R,a/x c , where realistic values of the di- 
mensionless parameter R are of the order k ~ 10~ 3 ; tak- 
ing a ~ (10 2 -j- 10 3 ) r c and x c ~ r c , we obtain k ~ 0.1 -j- 1 
which gives /3q ~ 3 -j- 13 according to Fig. [SJ 

V. FAR ZONE: MESO/MACROSCALE 
FRICTION 

At the mesoscopic scale, i.e. on distances r ^> A c , the 
substrate must be considered as deformable. Let us split 

I 



where we assumed that, for the sake of simplicity, the 
contacts are reborn with zero stretching, R(u) = 5{u), 
and 

r[X(r); r]= J di P(£) Qfc X(r); r] . (22) 

Equations (j2"Tl |2"21 should be completed with the elas- 
tic equation of motion for the sliding body (we assume 
isotropic slider) 

ii + ?/u = GiV 2 u + G 2 V(V-u), (23) 

where u(R) is the 3D displacement vector in the slider 
(R = {x,y, z}), 1] is the intrinsic damping in the slider, 
d = E/2(l + a)p = c\ and G 2 = G l /(l-2a)p = cf-cf, 
E, a and p are the Young modulus, Poisson ratio and 
mass density of the slider correspondingly, and c/ (c t ) is 
the longitudinal (transverse) sound speed. Equation (j2"3")l 
should be solved with corresponding boundary and initial 
conditions. In particular, at the interface (the bottom 
plane of the slider, where z = and {x, y} = r) we 
must have = X(r), u y = u z = 0, and the shear 
stress should equal FLY(r); r]/A 2 , where the friction force 
acting on the A c -block from the interface, 

F[X(r);r] = N x k c J duuQ[u;X(r);r], (24) 

should be obtained from the solution of Eq. (|2ip [here 
N x = (A c /a) 2 ]. 

Equations (|2"TH2"4"|) form the complete set of equations 
which describes evolution of the large scale tribological 



the frictional area into (rigid) blocks of size A c . In a 
general 3D model of the elastic slider, the nth A c -block 
is characterized by a coordinate X n , and its dynamics 
is described by the ME for the distribution functions 
Q n (u n ; X n ). A solution of these MEs gives the interface 
forces F n (X n ). Then, the transition from the discrete 
numbering of blocks to a continuum interface coordinate 
r is trivial: n -t r, Q n (u n ; X n ) -> Q[u; X(r);r], P n (u) -t 
P(u;r), T n {X n ) -► T[X(r);r], F n {X n ) F[X(r);r] 
(here r is a two-dimensional vector at the interface), and 
the master equation now takes the form: 



(21) 

I 

system; in a general case it has to be solved numerically. 
However, a qualitative picture may be obtained analyti- 
cally. The interface dynamics depends on whether or not 
the A c -blocks undergo the elastic instability, i.e., on the 
ratio of the stiffness of the A c -block K\ ps (2c 2 + 3c 2 )pA c 
[as follows from the discretized version of Eq. (l23l) ] and 
the effective critical stiffness parameter of the interface 
^Aeff = PoKx, where K* x - K Xs x c /Ax s and K Xs = 
N x k c . If the elastic instability does not appear, then 
a local perturbation at the interface relaxes, spreading 
over an area of size A s — the screening length considered 
below in Sec. IV Al In the opposite case, when the elas- 
tic instability does emerge (locally), in may propagate 
through the interface. Below in Sec. IV Bl we consider 
a simplified one-dimensional version, which allows us to 
get some analytical results and a rather simple simula- 
tion approach (such a model is also supported by the fact 
that the largest forces near the broken contact are just 
ahead/behind it according to Fig. [TJ. Recall that the 
interaction between the A c -blocks is weaker than in the 
short-range zone, it follows the law Sf cx r~ 3 which de- 
termines, e.g., the block-block interaction strength k x in 
Eq. (|2T)|) below (although the interaction is power-law, we 
may consider nearest neighbors only, because excitations 
at the interface, such as "kinks" introduced in Sec. IV Bl 
are localized excitations, and the role of long-range char- 
acter of the interaction reduces to modification of their 
parameters [24|). 



dQ[u;X{r);r] dQ[u; X(r);r] 

dXlr) du ^ U > ( " ^ = ^ ' r [ X( - r ) ; r l ' 
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A. Elastic screening length 



Let us assume that the slider is split in A c -blocks 
(rigid blocks) and consider the block-block interaction 
in a mean-field fashion (analogously to methods used in 
soft matter, see Refs. |25l - l27| ). Due to sliding of neigh- 
boring blocks, the forces acting on contacts in the nth 
A c -block get an addition shift. This effect may be ac- 
counted with the help of a substitution /„ — > f n + A/ n , 
A fn = Em/n/m x Prob(m -> broken) x U mn ps 
x c J2 m ^n fm F m Ti mn (recall that the sum is over the Ac- 
blocks here), or approximately 



1 + Xc X! T m(X m )U ri 



fn i 



(25) 



where T m (X m ) = J du P m {u) Q m {u; X m ) so that 
N\T rn (X m ) is the number of broken contacts in the mth 

I 



A c -block per its unit displacement, and 
n mn « N x n\ (A c /r mn ) 3 



(26) 



describes the dimensionless (i.e., normalized on f s ) elas- 
tic interaction between the A c -blocks separated by the 
distance r mn . In this way the force is given by / S II; the 
numerical constant k\ ~ Ra/X c depends on the substrate 
and interface parameters. 

Let us introduce the dimensionless variable e n = 
x c ^2 m , n T m (X m )I{ mn . The shift of forces in the nth. 
block due to broken contacts in the neighboring blocks 
may be accounted by a renormalization of the rate: 



Pn{u) -> P n [(1 + £»)«] 



(27) 



Indeed, when contacts in the neighboring blocks break, 
then the forces in the given block increase, e n > 0, and 
the contacts in the given block should start to break ear- 
lier, i.e., their threshold distribution effectively shifts to 
lower values. 

Making the transition from discrete sliding blocks to a 
continuum sliding interface, II mn — > II(r' — r) and e n — > 
e(r), we obtain a master equation of the form: 



dQ[u;X{r);r] dQ[u;X(r);r] 



dX(r) 



du 



P ([1 + e{r))u) Q[u- X(r);r] = 6(u) T[X(r);r] 



(28) 



where we again assumed that the contacts are reborn 
with zero stretchings, R(u) = d(u), 



e(r) = x c \- z I dVr[I(rV']n(r'-r) (29) 

-r|>Ac 



and 



T[X(r);r] = I dt; P ([1 + e(r)}0 Qfc; X(r);r] . (30) 

In the long- wave limit, when \de(r)/dr\ <C s(r)/\ c , 
we may assume that the interface is locally equili- 
brated, i.e., the distribution of forces on contacts is 
close to the steady-state solution of the master equation, 
Q[u; X(r);r] pa Q s (u;r), which depends parametrically 
on the coordinate r through the function e(r) entered 
into the expression for the rate P ([(1 + e(r)]u). The sta- 
tionary solution of the ME is known analytically [TTJ , and 
we may find the function (|3H| . T(r) = [1 + e(r)]/x c . To- 
gether with Eq. (j29|) this gives a self-consistent equation 
on the function e(r): 

e{r)=\- 2 f dV[(l + e(r')]II(r' - r) . (31) 

Taking into account the interaction of nearest neigh- 
boring A c -blocks only and expanding e(r) in Taylor series, 



we obtain the equation 



e(r) = n 



1 



l + e(r) + -\le"(r) 



(32) 



where IIo = ^II(A C ) = uN\K\ ~ vn\ c /a and v = 2 
4 is the number of nearest neighbors. Writing e(r) — 
e + Ae(r), where e = n /(l - IIo), Eq. (02) may be 
rewritten as 



X 2 s Ae"(r) = Ae(r) , 



(33) 



where A s = A c (eo/2) 1 / 2 is the characteristic screening 
length in the sliding interface. 

From the known analytical steady state solution of the 
ME [l2j], we may predict the dependence of screening 
length on temperature and sliding velocity. In particular, 
if T > 0, then X s oc w _1//2 — > oo as v — > in agreement 
with the results of Ref. El. 



B. Frictional crack as a solitary wave 

In the frictional interface, sliding begins at some weak 
place and then expands throughout the interface. Such 
a situation is close to the one known in fracture mechan- 
ics as the mode II crack, when the shear is applied along 
the fracture plane. In friction, a crack first opens, evolves 
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(propagates, grows, extends) during some "delay" time r, 
but then it either expands throughout the whole inter- 
face, or it will close because of the load. Below we con- 
sider the latter scenario, when one solid slips over another 
due to motion of the so-called self-healing crack [291432} 
— a wave or "bubble" of separation moving like a crease 
on rug [33j]. Our plan is to adopt ideas from fracture 
mechanics, adapt them to the friction problem and then 
reduce it to the Frenkel-Kontorova (FK) model [HJ in 
order to describe collective motion of contacts in the fric- 
tional interface. 

When one of the "collective contacts" (the A c -block) 
breaks, it may initiate a chain reaction, with contacts 
breaking domino-like one after another. This scenario 
may be described accurately by reducing the system of 
contacts to a Frenkel-Kontorova-like model. Recall that 
the FK model describes a chain of harmonically inter- 
acting atoms subjected to the external periodic potential 
V su h(x) of the substrate. If the atoms are additionally 
driven by an external force /, then the equations of mo- 
tion for the atomic coordinates u n take the form 

mu n + mriu n ~ g(u n+ i + u n -i - 2u n ) + F s ' ub (uj) = / , 

where m is the atomic mass, g is the strength of elas- 
tic interaction between the atoms, and rj is an effective 
damping coefficient which describes dissipation phenom- 
ena such as the excitation of phonons etc. in the sub- 
strate. The main advantage of using the FK model is 
that its dynamics is well documented [2J]. Mass trans- 
port along the chain is carried by kinks (antikinks) - 
local compressions (extensions) of otherwise commensu- 
rate structure. The kink is a well-defined topologically 
stable excitation (quasiparticle) characterized by an ef- 
fective mass rrik which depends on the kink velocity Vk, 
rnfe = m^l — w fe/c 2 ) -1 ^ 2 (the relativistic Lorentz con- 
traction of the kink width when its velocity approaches 
the sound speed c). Therefore, the maximal kink veloc- 
ity Wfcmax = c. In the discrete chain, kinks move in the 
so-called Peierls-Nabarro (PN) potential, whose ampli- 
tude is much lower than that of the primary potential 
V su b(x). Therefore, the kink motion is activated over 
these barriers, and its minimal velocity Vk m in is nonzero. 
The steady-state kink motion is determined by the en- 
ergy balance: the incoming energy (because of action 
of the external driving force /) should go to creation of 
new "surfaces" (determined by the amplitude of the sub- 
strate potential) plus excitation of phonons by the mov- 
ing kink (described by the phenomenological damping 
coefficient 77), so that Vk(f) = f /(rrikVi)- 

FK-ME model. Thus, let us consider a chain of Ac- 
contacts ("atoms" of mass m = pX^.), coupled harmoni- 
cally with an elastic constant g, driven externally through 
a spring of elastic constant K with the end moving with 
a velocity v. Using the discretized version of Eq. (f2"3"|). 
the elastic constants may be estimated as g 2X c pcf 
and K w \ c pc\. The A c -contacts are coupled "friction- 
ally" with the bottom substrate; the latter is described 
by the nonlinear force F s (u). The equation of motion of 
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FIG. 6. (color online): Color map of atomic velocities for a 
typical evolution of the chain of contacts. The nearest neigh- 
boring contacts interact elastically with the constant g = 25. 
The interaction with the substrate is modeled by the function 
F a (u) = fc c [tanh(ti) + 1.5e~" sin(3u)] with k c = 1 defined for 
< u < u c = 1 and periodically prolonged for other val- 
ues of u. All contacts are driven through the springs of the 
elastic constant K = 0.07, their ends moving with the veloc- 
ity v — 10 -4 . The motion is overdamped (m = 1, rj = 100). 
To initiate the breaking, two central contact interact with the 
substrate with smaller values of the elastic constant, k c — 0.5. 

the discrete chain is 

mu n +mr)u n -g(u n+ i+u n -i-2un)+F s (u n )+Ku n = f , 

where the driving force is given by f(t) — Kvt. The 
substrate force F s (u) is found from the solution of the 
ME for the rigid A c -block. A typical evolution of the 
chain is shown in Fig. [6l 

The general case may only be investigated numerically. 
Let us first consider a simplified case, when F s (u) has the 
sawtooth shape, i.e. it is defined as 

F s (u) — k c u for < u < u c (35) 

and periodically prolonged for other values of u. We 
assume that / is approximately constant during kink 
motion (otherwise, the kink will accelerate during its 
motion along the chain); this is correct if the change 
of the driving force A/ = Kv A< during kink motion 
through the chain, At — L/vk (L is the chain length 
and Vk is kink velocity), is much lower than k c u c , or 
K/k c < {v k /v)(u c /L). 

Let us define the function F(u) = F s (u)+Ku— f. The 
degenerate ground states of the chain are determined by 
the equation JF(u) = 0. Let the right-hand side (n — > 00) 
of the chain be unrelaxed, k c UR + Kur — /, or 

ur = f/{k c + K) , (36) 

while the left-hand side (n — > —00) already undergone 
relaxation, k c (iiL — u c ) + Kul = /, or 

u L = (f + k c u c )/(k c + K). (37) 
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Thus, the FK-like model of friction (the FK-ME 
model) is described by Eqs. (fMf and (j3"5)l with the bound- 
ary conditions given by Eqs. (|36|) and (f37|) . 

Continuum- limit approximation. Let the system be 
overdamped (ii = 0); later on we shall remove this restric- 
tion. In the continuum-limit approximation, n — > x = na 
(a = 1), the motion equation takes the form 

mrjut - a 2 gu xx + F(u) = 0, T(u)\ x ^ ±oa = . (38) 

We look for a solution in the form of a wave of sta- 
tionary profile (the solitary wave), u(x,t) = u(x — Ufet), 
so that Ut — —Vku' and u xx = u" . In this case Eq. 
takes the form 



mrjVku' + a 2 gu" = Fiu) , 



(39) 



which may be solved analytically by standard meth- 
ods (HI- 

A solution of Eq. (|39l) with these boundary conditions 
exists only for a certain value of the kink velocity Vk, 
defined by the equation 



(m V v k ) 2 = ga 2 (k c + K){2 - /3) 2 /(/3 - 1) 



(40) 



where /3 — fc c /(fc* — K) and fe* = f/u c . The solitary- 
wave solution exists for forces f min < f < / max only. 
The minimal force which supports the kink motion - 
the Griffith threshold — is given by 



fn 



2 kc 



The maximal force, for which a kink may exist, is given 
by 



/max — (kc + Kj U c ; 



(42) 



at higher forces, the barriers of F s (u) are degraded, the 
stationary ground states disappear, and the whole chain 
must switch to the sliding state. 

From Eq. (l40l) we can find the kink velocity as a func- 
tion of the driving force. At low velocities 

Vk » (/ - /min)/"lfe»7 ) (43) 

where we introduced the effective kink (crack) mass 



mk=m /vM, ( 1+ f ' 



while at / — > / max the velocity tends to infinity, 



mrjVk 



I gk c (k c + K)a 2 u c 

(/max /) 



(44) 



(45) 



The latter limit should be corrected by taking into ac- 
count inertia effects. The term mil in Eq. (|34[) gives 



mv 2 u" for the solitary-wave solution, so it can be in- 
corporated if we substitute in the above equations g — > 
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FIG. 7. (color online): Evolution of the chain of N = 3000 
contacts. The nearest neighboring contacts interact elasti- 
cally with the constant g = 25, the interaction with the sub- 
strate is modeled by the sawtooth function ([35} with fc c = 1 
and u c = 1. All contacts are driven through the springs of the 
elastic constant K = 0.07, their ends moving with the veloc- 
ity v — 10~ 4 . The motion is overdamped (m = 1, rj = 100). 
To initiate the breaking, two central contacts interact with 
the substrate with smaller spring constants, k' c = 0.5. When 
the kinks motion begins, the elastic constants of the central 
contacts restore their values to k c = 1, and the driving ve- 



(41) 

locity changes its sign, v 



v b = -2 x 10" 



(a) shows the 



kinks centers (defined as places where the atomic velocity is 
maximal), (b) shows the driving force f(t), (c) shows the av- 
erage chain velocity = iV" 1 J^. in, and (d) demonstrates 
oscillation of the velocity due to PN barriers. 



9cS = 5(1 — v t/ c o)i where cq = (ga 2 /m) 1 / 2 is the sound 
speed along the chain. The high- velocity limit now takes 
the form 




(46) 



Simulations. The continuum-limit approximate is ac- 
curate for the case of strong interaction between the con- 
tacts, g ^> 1; in the opposite limit one has to resort to 
computer simulation. We solved Eq. (j34|) by the Runge- 
Kutta method. As the initial state, we took the chain 
of length N (typically N = 3 x 10 3 or 3 x 10 4 ) with 
periodic boundary conditions and all contacts relaxed, 
but the threshold breaking value for two central contacts 
was taken lower than for the other contacts. Then the 
driving force increases because of stage motion, two cen- 
tral contacts break first and initiate two solitary waves 
of subsequent contact breaking which propagate in the 
opposite directions through the chain. The value k' c of 
the lower threshold of the central contacts determines the 
driving force and therefore the kink velocity; the lower 
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FIG. 8. (color online): Kink velocity versus the driving force 
for (a) g = 5 (v b = -4 x 1(T 5 ) and (b) g = 25 (v b = -2 x 
10~ 4 ); N = 3 x 10 4 , other parameters as in Fig. [7J Blue 
solid and red dashed lines correspond to Eqs. (|40l) and (|43|) . 
correspondingly. 



this threshold, the lower the threshold force for the mo- 



tion to start 3J|. As soon as the kink motion is initi- 
ated, the k c - values of the central contacts are restored 
to the same value as for other contacts (otherwise these 
contacts will act as a source of creation of new pairs of 
kinks), and we begin to move the stage in the opposite 
direction, v > — > V(, < 0, so that the driving force lin- 
early decreases with time (see Fig.[7b>), the average chain 
velocity («,) = iV -1 in decreases as well (Fig. [7b) un- 
til the motion stops (Fig. UK)- Also, such an algorithm 
allows us to find the dependence of the kink velocity de- 
termined as 



Vk 



= n k 1 N({u l )-v), 



(47) 



where — 2 is the number of moving kinks in the chain 
and v = ul.r — VbK / (k c + K) is the background velocity, 
on the driving force /. These dependences are presented 
in Fig. [5] they agree well with that predicted by Eqs. P0|) 
and g3]). 

Contrary to the continuum-limit approximation, in the 
discrete chain of contacts the kink oscillates during mo- 
tion (see Fig. [7bl) — the well-known discreteness effect of 
the FK model due to existence of the PN barriers /pn- 
The stronger the elastic interaction between the contacts, 
the larger the kink "width" and the smaller the kink oscil- 
lations (compare Figs.[5£i and |SJd). The amplitude of os- 
cillations also depends on the shape of the "substrate po- 
tential" [24| — it is larger for a sawtooth potential F s (u), 
but smaller for a smoother shapes. Recall that the Ac- 
contacts are characterized by a smooth dependence F s (u) 
as follows from the master equation. The PN oscillations 



determine the lowest average kink velocity. Therefore, 
the lowest velocity allowed for the frictional crack propa- 
gation, Vk mm, is determined by the parameters g and A c 
— the larger are g and A c , the smaller is vt min- 

Discussion. The FK-ME model used here is rather 
close to the well-known ID Burridge-Knopoff (BK) 
model of earthquakes with a velocity-weakening friction 
law (35[. The difference is in the interface force F s (u): we 
use the function derived from the ME-EQ model (with 
well-defined parameters which may be extracted from ex- 
periments or calculated from first principles), whereas the 
BK model adopts a phenomenological velocity-dependent 
function for F s . Nevertheless, the qualitative behavior 
of the two models is similar, the BK model also exhibits 
solitary-wave dynamics as was demonstrated numerically 
in Ref. [36j . In our case, however, by reducing the model 
to the FK-ME one, we can describe the solitary waves 
analytically and rigorously. 

In the simulation we started from the well-defined ini- 
tial configuration, when all contacts are relaxed except 
the one or two where kink's motion is initiated. If one 
starts from a random initial configuration, we expect that 
kinks will emerge at random places, so that several kinks 
may propagate through the system simultaneously, as 
was observed in simulation of the BK model [36| . 

Also we assumed that all A c -contacts are characterized 
by the same F s (u) dependence and thus have the same 
threshold values i^ii- This is correct if the number of orig- 
inal contacts within a single A c -contact, N\ = (A c /a c ) , is 
infinite. Otherwise, different A c -contacts will have differ- 
ent threshold values F n however the distribution of their 
thresholds is narrower that the distribution of thresholds 
of single asperities by a factor ^/Nx. A narrow distribu- 
tion of thresholds will nevertheless have a qualitative ef- 
fect because rupture fronts may stop when they meet Ac- 
contacts with a threshold above the driving force. When 
the interface is disordered, the avalanches will have fi- 
nite lengths and may become short for forces near /i n j, 
for which the rupture fronts propagate at the minimal 
velocity. 

Our approach may also incorporate the existence of 
disorder and defects always present in real materials. On 
the one hand, defects may nucleate kinks (cracks); on the 
other hand, the kink propagation may be slowed down 
up to its complete arrest due to pinning by the defects. 
For example, the slowing down of the ID crack prop- 
agating through a 2D system with quenched randomly 
distributed defects was considered in Ref. |37| . 

Thus, reducing the EQ-ME model of friction to the 
FK-ME one, we described avalanche- like dynamics of the 
frictional interface — the solitary wave of contacts break- 
ing. If the force F s (u) has a sawtooth shape, then the 
interface dynamics may be described analytically; other- 
wise one has to use numerics. The analogy with the FK 
model may be extended even further: 
• The driven FK model exhibits hysteresis when the force 
increases and then decreases [U, [38| . The same effect 
was observed in the large-scale crack simulation [39( , thus 
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could be observed in the frictional interface too. 

• Effects of nonzero temperature may be considered. One 
may predict that at T > the sliding kinks will expe- 
rience an additional damping, while the immobile (e.g., 
arrested) kinks will slowly move (creep) due to thermally 
activated jumps. 

• As shown in Refs. [24], E3|, a fast driven kink begins 
to oscillate due to excitation of its shape mode, and 
then, with the further increase of driving, the kink is 
destroyed. This effect is similar to what is observed in 
fracture mechanics, where cracks begin to oscillate and 
then branch [4l| . 

• If the interaction between the atoms is nonlinear and 
stiff enough, the FK model admits the existence of super- 
sonic kinks [42| which are similar to solitons of the Toda 
chain. It would be interesting to study if similar waves 
may appear in the frictional interface, as was predicted 
in crack propagation (43|. 

• One may suppose that the damping coefficient r\ in the 
equation of motion (f3"4"|) depends on the kink velocity, 
Tj(v). In fracture mechanics, this coefficient defines the 
rate at which the energy is removed from the crack edge, 
thus it plays a crucial role. 

• A large number of works is devoted to different gener- 
alizations of the FK model to 2D system (e.g., see [Hj]). 
For example, if kinks attract one another in the y (trans- 
verse) direction, they unite into a line (dislocation) which 
moves as a whole (or due to secondary kinks). 

VI. CONCLUSION 

We discussed the crucial role in sliding friction of 
the elastic interaction between the contacts at the inho- 
mogeneous frictional interface and proposed various ap- 
proaches to treat this problem from different viewpoints. 
The interaction produces a characteristic elastic correla- 
tion length X c — a 2 E/k c . At distance r < X c the slider 



may be considered as a rigid body but with a strong 
contacts' interaction, which leads to shrinking of the ef- 
fective contact breaking threshold distribution and an 
enhanced possibility for a mechanical elastic instability 
to appear, which is conducive to stick slip. At large dis- 
tances r > A c , the contact-contact interaction leads to 
screening of local perturbations in the interface, or to 
appearance of collective modes (frictional cracks) propa- 
gating as solitary waves. 

In our work we assumed that the external stress (the 
driving force) is uniform across the system. In a gen- 
eral case, however, stress is nonuniform and may more- 
over change with (adjust itself to) interface dynamics, so 
that the problem should be considered self-consistcntly. 
For given boundary conditions, determined by the exper- 
imental setup, one should calculate the stress field, e.g., 
by finite element technique, which provides the driving 
force f(r) in the FK-ME model. The latter defines the 
displacement field at the interface through the solution 
of the FK-ME master equations. The displacement field 
in turn is to be used as the boundary condition for the 
elastic-theory equations at the frictional interface (from 
other sides of the slider, the boundary conditions should 
correspond to a given experimental setup). 
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